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1.  Introduction 


Duct  flows  occur  in  a  variety  of  engineering  applications.  They  are  present 
in  aircraft  and  automobile  engines,  heat  exchangers,  steam  generators,  and  piping 
systems.  Often,  the  ducts  have  a  noncircular  cross  section,  which  undergoes  a 
change  of  shape  and  size  along  the  length  of  the  duct.  The  purpose  of  the  present 
research  project  is  to  develop  a  calculation  procedure  for  the  prediction  of  flow 
in  such  ducts. 

Many  duct  flows  are  characterized  by  the  absence  of  reverse  flow  or  separa¬ 
tion  and  by  a  nearly  uniform  pressure  over  any  cross  section.  These  conditions 
are  found  when  the  change  in  the  size  and  shape  of  the  cross  section  is  gradual 
rather  than  sudden.  Such  flows  can  be  treated  as  fully  parabolic,  i.e.,  they  can 
be  predicted  by  a  marching  procedure  starting  at  the  inlet  plane  and  proceeding 
to  successive  cross-sectional  planes  downstream.  The  procedure  described  in  this 
report  is  of  the  fully  parabolic  variety. 

If  significant  pressure  variations  arise  over  a  cross  section,  they  can  be 
incorporated  in  a  procedure  known  as  partially  parabolic.  In  such  a  procedure, 
the  duct  length  must  be  swept  many  times  by  the  marching  calculation  until  the 
three-dimensional  pressure  variations  are  properly  established  and  used  in  the  flow 
calculation.  If  a  reverse  flow  Is  present  in  the  duct,  a  fully  elliptic  procedure 
is  needed  for  its  calculation. 

The  fully  parabolic  procedure  employed  in  this  report  calculates  the  flow  over 
one  cross  section  at  a  time.  Thus,  the  computational  task  is  almost  the  same  as 
solving  a  two-dimensional  problem  over  a  domain  of  arbitrary  shape.  The  method  of 
Ballga  and  Patankar  [1,  2]  was  chosen  to  be  the  starting  point  for  this  purpose. 

In  the  course  of  adapting  that  method,  an  Improvement  in  the  treatment  of  pressure 
was  worked  out.  The  two-dimensional  calculation  scheme  along  with  the  proposed 
improvement  is  described  in  Section  2  of  this  report.  The  method  has  been  tested 


on  a  large  number  of  problems  to  evaluate  Its  correctness,  accuracy,  convergence, 
and  other  behavior.  Results  of  some  of  the  test  problems  are  presented  In  Section 
3. 

The  extension  of  the  two-dimensional  method  to  three-dimensional  duct  flows 
is  described  in  Section  4.  Section  5  is  devoted  to  some  applications  of  the 
method.  They  include  the  developing  flew  in  some  constant-area  ducts;  but  of 
particular  Importance  are  the  calculations  for  two  ducts  of  varying  cross-sectional 
area  and  shape. 
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2.  Method  for  two-dimensional  flows 


Uhen  a  three-dimensional  duct  flow  is  treated  as  fully  parabolic  in  nature, 
its  calculation  is  performed  by  a  marching  procedure  in  the  main  flow  direction. 

At  each  marching  station,  a  computationally  two-dimensional  problem  is  solved  over 
the  duct  cross  section.  Thus,  a  procedure  for  two-dimensional  flows  forms  an 
essential  building  block  of  the  method  for  three-dimensional  parabolic  flows.  The 
two-dimensional  procedure  employed  in  the  present  work  is  described  in  this  section. 

It  was  Intended  to  employ  the  two-dimensional  finite-element  procedure  of 
Ballga  and  Patankar  [1,  2]  for  this  building  block.  However,  the  treatment  of  the 
pressure  field  in  [2]  is  somewhat  Inconvenient;  therefore,  an  improved  treatment  was 
worked  out.  It  is  this  improved  procedure  that  is  described  here. 

2.1  Governing  differential  equation 


As  described  in  [3] ,  the  fluid  flew  and  heat  transfer  phenomena  in  steady  state 
are  governed  by  a  general  differential  equation  of  the  form: 

■£r}  <«v>  - <r  +  s  <2-i) 

where  <J>  is  the  general  dependent  variable,  pu^  denotes  the  mass  flow  rate  causing 
the  convective  transport,  T  is  the  diffusion  coefficient,  and  S  represents  the  source 
term.  The  variable  $  can  stand  for  enthalpy  (or  temperature),  concentration,  turbu¬ 
lence  energy,  or  a  velocity  component.  Thus,  a  solution  procedure  for  any  convective 
phenomenon  should  provide  a  mechanism  for  solving  Gq.  (2.1). 

2.2  Domain  discretization 

In  the  calculation  procedure  described  in  [1,  2],  the  physical  domain  of  inter¬ 
est  is  discretized  into  triangular  elements  as  shown  in  Fig.  2.1.  The  vertices  of 
the  triangles  are  called  nodes  or  grid  points.  The  aim  of  the  calculation  scheme  is 
to  obtain  the  values  of  <p  at  these  nodes. 


Fig.  2.1  Domain  discretization 


The  algebraic  equations  for  the  $  values  at  the  nodes  are  obtained  by  integrating 
Eq.  (2.1)  over  the  polygonal  control  volume  associated  with  each  node.  The  construc¬ 
tion  of  these  control  volumes  is  also  shown  in  Fig.  2.1.  Each  triangular  element  is 
split  into  three  parts  by  joining  its  centroid  to  the  centers  of  the  three  sides.  The 
control  volume  around  a  typical  node  P  is  composed  of  the  parts  of  all  the  neighboring 
triangles  that  contain  the  point  P. 

2.3  Discretization  equation 

The  discretization  form  of  Eq.  (2.1)  is  obtained  by  integrating  the  differential 
equation  over  a  control  volume.  Such  integral  equation  consists  of  the  total  source 
term  S  over  the  control  volume  and  the  net  convecton  and  diffusion  flux  across  the. 
faces  of  the  control  volume.  As  seen  from  Fig.  2.1,  the  control-volume  faces  are 
composed  of  the  dashed  lines  that  lie  within  each  triangle.  To  express  the  total 
flux  (i.e. ,  the  convection  flux  plus  the  diffusion  flux)  across  such  a  dashed  line 
in  terms  of  the  nodal  values  of  a  shape  function  is  required. 

As  explained  in  [1],  an  appropriate  shape  function  for  convection-diffusion 
problems  is  one  that  takes  account  of  the  resultant  flow  direction.  The  shape 


function  suggested  in  [1]  and  implemented  in  this  work  is 


<t>  -  A  +  B  exp(pUX/D  +  CY  (2.2) 

where,  as  shown  in  Fig.  2.2,  the  coordinate  X  is  aligned  in  the  resultant  flow  dir¬ 
ection,  while  Y  is  taken  as  normal  to  it.  The  quantity  U  represents  the  magnitude 
of  the  velocity  vector.  The  constants  A,  B,  and  C  in  Eq.  (2.2)  are  expressible  in 
terms  of  <j>^,  >  <t>?  at  the  three  nodes  of  the  triangle. 


Fig.  2.2  Coordinate  axes  for  shape  function 
The  use  of  Eq.  (2.2)  in  the  integral  equation  for  the  control  volume  leads  to 
the  algebraic  equation 


*p 


Ea  ,  <p  .  +  b 
nb  Tnb 


(2.3) 


where  the  subscript  nb  denotes  a  neighbor  node  of  P  and  the  summation  is  taken 
over  all  the  neighbors  of  P.  The  coefficients  ap  and  a^  represent  the  convection- 
diffusion  Influence,  while  b  contains  the  integral  of  the  source  term  S  over  the 
control  volume. 


2.4  Solution  of  he  discre'  .zatlon  equations 
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of  the  proposed  method,  the  domain  discretization  is  performed  by  the  use  of  a 
line  structure.  Although  this  places  seme  restriction  on  the  shape  of  the  domain 
or  on  the  local  grid  fineness,  the  technique  provides  considerable  convenience  in 
node  numbering  and  in  the  solution  of  equations.  With  the  nodes  arranged  on  a 
line  pattern,  the  line-by-line  method  [3]  becomes  the  natural  choice  for  solving 
the  algebraic  equations.  It  is  this  method  that  is  employed  in  the  present  work  . 

2.5  Calculation  of  fluid  flow 

Although  the  general  dependent  variable  <t>  can  stand  for  the  two  velocity  com¬ 
ponents  u  and  v  in  a  two-dimensional  flow  and  thus  a  calculation  procedure  for  fluid 
flow  is  contained  within  the  general  procedure  for  solving  Eq.  (2.1),  there  is  an 
important  complication.  The  source  term  S  contains  the  pressure  gradient  9p/9x  or 
9p/9y  when  <f>  stands  for  u  or  v.  The  pressure  field  p  represents  an  extra  degree  of 
freedom  and  it  is  indirectly  constrained  by  the  requirement  that  the  velocity  field 
must  satisfy  the  continuity  equation 

9  (Pu.)  -  0  (2.4) 

3xi 

The  combined  treatment  of  the  momentum  and  continuity  equations  forms  the  special 
features  of  a  fluid-flow  calculation. 

If  the  velocity  components  and  pressure  are  stored  at  the  same  grid  points, 
then  the  required  interpolations  of  pressure  and  velocity  lead  to  discretization 
equations  that  admit  unrealistic  zig-zag  velocity  fields  and  checkerboard  pressure 
fields.  This  is  explained  at  length  in  [3].  The  remedy  adopted  in  finite-difference 
methods  is  a  staggered  grid,  in  which  the  velocity  components  are  stored  at  locations 
that  are  displaced  (or  staggered)  relative  to  the  pressure  nodes.  The  same  remedy 
is  not  available  in  finite  element  methods,  since  the  lines  joining  adjacent  nodes 
do  not  lie  along  the  coordinate  directions.  The  common  remedy  is  a  mixed-inter¬ 
polation  (or  unequal-order)  scheme,  in  which  pressure  is  stored  at  fewer  locations 


than  the  velocity  components.  This  practice  is  employed  in  (2]  and  also  by  earlier 
workers . 

Although  the  unequal-order  procedure  is  adequate,  it  is  not  entirely  satis¬ 
factory.  It  leads  to  complications  in  numbering  sequences,  control-volume  defini¬ 
tions,  etc.,  since  the  pressure  grid  and  the  velocity  grid  must  be  handled  separa¬ 
tely.  Moreover,  the  coarser  grid  used  for  pressure  would  limit  the  accuracy  of  the 
solution.  If  the  pressure  grid  is  made  sufficiently  fine,  the  associated  velocity 
grid  becomes  excessively  fine  and  thus  wasteful.  In  the  interests  of  convenience 
and  efficiency,  an  equal-order  procedure  is  thus  sought. 

2.6  Discretized  momentum  equations 

Since  the  momentum  equations  are  special  cases  of  the  general  differential 
equation  (2.1),  they  can  be  cast  into  the  discretization  form  similar  to  Eq.  (2.3). 
Thus 


a 

u  =  E 

a  , 

U  , 

+ 

b 

P 

P 

nb 

nb 

a 

v  -  £ 

a  , 

V  , 

+ 

b 

P 

P 

nb 

nb 

U  V 

where  b  and  b  are  the  integrals  (over  of  the  control  volume) 
(Su  -  9p/  9x)  and  (Sv  -  9p/9y) 


(2.5) 

(2.6) 

of  the  source  terms 


tespectively. 

For  further  rearrangement  of  the  equations,  pseudovelocities  u  and  v  are 
defined  by 

%  ~  (E  anb  unb)/aP  (2 


and 


*P  '  (Eanb  v„b> /*P 

Then  the  momentum  equations  (2.5)  and  (2.6)  can  be  rewritten  as 

Up  ■  Up  +  du  (SU  -9p/9x)c  v 


(2.8) 


(2.9) 


(2.10) 


vp  *  vp  +  dv  (sv  -  3p/9y>C  Vi 

where  the  subscript  c.v.  indicates  an  average  value  over  the  control  volume.  The 

U  V 

multipliers  d  and  d  can  be  seen  to  be 

du  -  dV  =  AV/ap  (2.11) 

where  AV  is  the  volume  of  the  control  volume. 

No  novelty  has  yet  been  introduced.  If,  in  the  continuity  equation,  a  linear 
interpolation  of  u  and  v  is  used,  the  well-known  difficulties  associated  with 
storing  all  variables  at  every  node  would  result.  The  equal-order  procedure  here 
is  developed  by  proposing  a  different  interpolation  technique  for  the  continuity 
equation.  This  technique  is  described  next. 

2.7  Treatment  of  the  continuity  equation 

The  discretized  form  of  Eq.  (2.4)  would  contain  mass  flow  rates  across  the 
boundaries  of  a  typical  control  volume  shown  in  Fig.  2.1.  These  flow  rates  can  be 
compiled  by  calculating  the  flow  crossing  the  three  dashed  line  segments  in  every 
element.  A  typical  element  is  shown  in  Fig.  2.3  . 


3 


Fig.  2.3  Control-volume  faces  within  an  element 
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Although  it  is  possible  to  calculate  u  and  v  at  any  point  within  the  element 

by  a  linear  interpolation  between  the  nodal  values  of  u  and  v,  the  practice  proposed 

here  makes  the  velocities  appearing  in  the  continuity  equation  directly  responsive 

to  the  pressure  gradient  over  the  element.  Thus,  for  the  purpose  of  the  continuity 

equation,  the  velocity  field  within  an  element  is  calculated  from  u  and  v,  where 

u  -  u  +  dU  [SU  -  3p/3x]e  (2.12) 

v  -  v  +  dv  [SV  -  3p/3y]fi  (2.13) 

where  the  subscript  e  indicates  the  element  under  consideration.  The  source  terms 
u  v 

S  and  S  are  considered  to  be  uniform  over  an  element.  Further,  if  the  pressure  p 
is  described  by  a  linear  shape  function  over  the  element,  the  gradients  3p/3x  and 
3p/3y  also  become  uniform.  The  calculation  of  u,  v,  d^and  dV  at  the  nodes  has 
been  outlined  in  Eqs.  (2.7),  (2.8),  and  (2.11).  For  the  calculation  of  u  and  v 
from  Eqs.  (2. 12) -(2. 13) ,  the  values  of  u,  v,  dU,  and  dV  are  considered  to  vary 
linearly  over  the  element.  As  a  result,  the  u  and  v  values  would  have  a  linear 
distribution. 

That  the  u  and  v  field  as  given  by  Eqs.  (2. 12)-(2.13)  is  "driven"  directly  by 
the  nodal  pressures  in  a  given  element  is  the  essential  feature  of  this  formulation. 
This  eliminates  the  possibility  of  checkerboard  pressure  fields. 

To  obtain  a  discretization  equation  for  pressure,  first  a  discretized  form  of 
the  continuity  equation  is  written  by  integrating  it  over  the  typical  control 
volume  shown  in  Fig.  2.1.  The  mass  flow  rates  across  the  faces  of  the  control 
volume  are  expressed  in  terms  of  u  and  v  defined  in  Eqs.  (2.12)  -  (2.13).  The 
pressure  gradients  3p/3x  and  3p/3y  are  written  in  terms  of  the  nodal  pressure. 

The  result  is  the  following  discretization  equation  for  p,  which  has  the  same  appear¬ 
ance  as  Eq.  (2.3) 

»P  Pp  ■  1  anb  pnb  +  b  <2'U) 


-  1 
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where  P  denotes  a  typical  node  shown  In  Fig.  2.1,  nb  stands  for  a  neighbor  node 

of  P,  and  the  summation  Is  taken  over  all  the  neighbors.  The  coefficients  a^  and 

a  ,  arise  from  the  multipliers  dU  and  dv  in  Eqs.  (2.12)-(2.13) .  The  term  b  in 
nb 

Eq.  (2.14)  contains  u,  v,  SU,  and  Sv.  Since  the  structure  of  Eq.  (2.14)  is  identical 
to  that  of  Eq.  (2.3),  the  line-by-line  solution  method  mentioned  in  Section  2.4  is 
applicable  to  Eq.  (2.14)  as  well. 

An  alternative  form  of  Eq.  (2.14)  can  be  derived  by  defining  a  pressure  cor¬ 
rection  p'  such  that 

p  «p*  +  p*  (2.15) 

where  p*  is  the  current  estimate  or  guess  for  pressure.  If  the  momentum  equations 
(2.5)  -  (2.6)  are  solved  by  substituting  the  guessed  pressure  field  p*,  the  resulting 
velocities  can  be  denoted  by  u*  and  v*  .  Consequently,  a  velocity  field  u*  and  v* 
can  be  obtained  from  Eqs.  (2.12)  -  (2.13)  by  substituting  p*  for  p  and  by  calculating 
u  and  v  from  the  tf*and  v*values.  Then  the  counterpart  of  Eq.  (2.14)  can  be  written 
as 

*P  PP  "  Z  anb  pnb  +  b  (2,16) 


where  ap  and  a^  are  Identical  to  the  corresponding  coefficients  in  Eq.  (2.14)  and  b 

*v  ^  JU 

results  from  the  substitution  of  the  u  and  v  field  in  the  continuity  equation. 


2.8  The  overall  solution  procedure 

Since  all  the  ingredients  of  the  calculation  procedure  have  been  outlined,  the 
overall  procedure  can  now  be  described  in  terms  of  the  various  steps  in  the  calcu¬ 
lation  sequence. 


(1)  Start  with  initial  guesses  for  u,  v,  u,  v,  and  other  relevant  variables. 

(2)  Calculate  the  coefficients  in  the  momentum  equations  for  u  and  v. 

(3)  Hence  obtain  u,  v,  du,  and  dVfrom  Eqs.  (2.7),  (2.8),  and  (2.11). 

(4)  Set  up  Eq.  (2.14)  and  solve  it  to  get  a  pressure  field. 

(5)  Treating  this  as  p*,  solve  the  momentum  equations  to  obtain  u*  and  v* 


(6)  Hence  set  up  the  pressure-correction  equation  (2.16)  and  solve  it. 

Use  the  resulting  p'  values  to  correct  u* ,  v*,  u*  and  v*  .  At  this  stage,  the 
corrected  field  u,  v  would  satisfy  continuity. 

(7)  Solve  Eq.  (2.3)  for  other  relevant  dependent  variables  such  as  tempera 
ture,  concentration,  turbulence  parameters,  etc. 

(8)  Return  to  step  2  and  repeat  until  convergence. 
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3.  Application  of  the  method  to  some  two-dimensional  problems 

In  this  section,  the  proposed  method  is  applied  to  three  test  problems. 

3.1  Flow  between  two  concentric  rotating  cylinders 

The  flow  of  an  incompressible  fluid  between  two  concentric  rotating  cylinders 
is  shown  schematically  in  Fig.  3.1.  The  inner  cylinder  is  supposed  to  be  at  rest 
while  the  outer  cylinder  rotates  with  angular  velocity  w.  In  the  (r-0)  polar 
coordinates,  this  is  a  one-dimensional  problem  with  the  radial  component  of  velo¬ 
city  ur  being  identically  zero.  However,  in  the  Cartesian  coordinates  (x-y), 
the  problem  is  fully  two-dimensional,  with  non-vanishing  velocity  components  u  and 
v  along  the  x  and  y  coordinates  respectively.  The  aim  here  is  to  solve  this  two* 
dimensional  problem  in  the  x-y  coordinates  over  a  square  region  R  shown  shaded 
in  Fig.  3.1. 

The  exact  solution 

Let  Uq,  ur  and  p  be,  respectively,  the  circumferential  velocity,  the  radial 
velocity  and  pressure.  Also  let  u  and  v  be  the  velocity  components  along  the 
x  and  y  directions.  At  this  point,  the  following  dimensionless  variables  are 
introduced  . 

*  U9 


u9  “  2r]U 

(3.1) 

*  .  Ur 

ur  “  2rju> 

(3.2) 

*  (P'Pq* 

P  "  2 

(3.3) 

p(2r  j  w) 

*  u 

(3.4) 

u  mw^r 

*  V 

(3.5) 

v  "T2i vT 

* 

r  * 


r 

rl 


(3.6) 
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* 

x 


(3.7) 


(3.8) 


where,  p  is  the  fluid  density,  and  pq  is  the  reference  pressure  on  inner  cylinder 
(r*r^).  For  a  radius  ratio  of  2  between  the  outer  and  inner  cylinders,  as  shown 
in  Fig.  3.1,  the  exact  solution  to  the  problem  is  given  by 


and 


r" 


0 


* 

P 


1,  *2 
9(r 


from  which  it  follows  that 


and 


* 


-¥‘ 


(3.9) 

(3.10) 


(3.11) 


(3.12) 


(3.13) 


Domain  discretization 

Figure  3.2  shows  the  discretization  of  the  computational  domain  into  quadri¬ 
laterals  and  further  subdivision  of  these  quadrilaterals  into  triangular  elements. 
As  can  be  seen,  the  uniform  N  x  N  grid  has  a  spacing  Ax*  ”  Ay*  m  6,  where. 


6 


Results  were  obtained  for  four  different  grid  sizes  corresponding  to  N  « 


(3.14) 
9,  13,  17 


and  21. 


Results 


(i)  Table  3-l(A)  shows  the  variation  of  percentage  error  in  u  at  the  mid 
point  of  the  computational  domain  with  Reynolds  number  and  grid  size.  T>.e  percent¬ 
age  error  is  defined  as: 

iU* 


-u 


percentage  error  in  u*  ■  100 


exact  computed 


u* 


exact 


(3.15) 


where  u* 


is  the  exact  value  of  u*  given  by  Eq.  (3.12)  and  u* 

computed 


exact 

is  the  numerically  confuted  value 


As  can  be  seen  in  Table  3-l(A) ,  the  percentage  error  increases  with  Reynolds 
number  at  any  given  grid  size.  For  any  given  Reynolds  number,  the  error  decreases 
continuously  as  the  grid  is  refined.  Table  3-1 (B)  shows  the  results  of  unequal 
order  method  of  Baliga  and  Patankar  [2]  for  a  21  x  21  grid.  Comparison  of  these 
with  the  present  results  for  a  21  x  21  grid  (Table  3-l(A))  shows  that  the  proposed 
method  gives  lower  error  or  better  results. 

(ii)  Let  e  be  the  average  of  the  percentage  error  defined  in  Eq.  (3.15). 
Further,  let 

e  -  a  6P  (3.16) 


where  a  is  a  constant  and  p  is  the  so-called  'order  of  convergence'.  If  z ^  and 


^2  are  the  values  of  E  for  grid  size  6^  and  62  respectively,  p  can  be  determined 
from 

1n(ej)-!n(e2) 


(3.17) 


P  “  177(6,  )-ln<5i) 

Table  3-2 (A)  shows  the  variation  of  E  with  grid  size  and  Reynolds  number  for 


the  equal  order  method.  Also  included  in  this  table  are  the  values  of  order  of  con¬ 
vergence  p  obtained  using  the  results  of  17  x  17  and  21  x  21  grid.  As  can  be  seen, 
the  percentage  error  over  the  domain,  e  increases  with  Reynolds  number  for  any  fixed 
grid  size.  Also,  at  any  fixed  Reynolds  number  e  decreases  with  increasing  grid 
size  confirming  that  the  computed  solution  converges  to  the  exact  solution  as  the 
grid  Is  refined.  From  the  tabulated  values  of  p,  it  may  be  inferred  that  the  order 


Re 

Grid 

1 

10 

100 

1000 

9  *  9 

.0213 

.  149 

.158 

<*> 

00 

13  x  13 

.0107 

.067 

.268 

.342 

17  x  17 

.00618 

.0378 

.238 

.206 

21  x  21 

.00404 

.0241 

.193 

.149 

TABLE  3-1 (B) 

Percentage  error  in  u*  at  the  center  of  domain  -  unequal  order  method 


Re 

Grid 

1 

10 

100 

1000 

21  x  21 

.02 

.034 

.264 

CO 

CM 
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TABLE  3-2 (A) 

Average  percentage  error  in  u*  over  the  domain  -  equal  order  method 


Re 

Grid;  8 

1 

10 

100 

1000 

9  x  9;  .l25/./f 

.0609 

.0952 

.337 

.637 

13  x  13;  .0833/vT 

.0251 

.0442 

.222 

.345 

17  x  17;  .0625/ 

.0135 

.025 

.157 

.226 

21  x  21;  .05/ 

.00846 

.0162 

.114 

.166 

Order  of  converg¬ 
ence  p,  using 
results  of  1 7X 1 7 
21x21  grids 

2.094 

1-944 

1.434 

1.383 

TABLE  3*2 (B) 

Average  percentage  error  in  u*  over  the  domain  -  unequal  order  method 


Re 

Grid; 8 

1 

10 

100 

1000 

17  x  17;  .0625//T 

.050 

.093 

.545 

.803 

21  x  21;  .05//T 

.034 

.060 

.407 

.594 

Order  of  con¬ 
vergence  p 

1.728 

1.964 

1.308 

1.351 

of  convergence  of  the  method  is  between  1  and  2. 


Table  3-2 (B)  shows  the  corresponding  results  for  the  unequal  order  method  of 
Baliga  and  Patankar  [2]. 

(iii)  Fig.  3.3(a)  shows  a  comparison  between  the  exact  and  computed  solution 
for  pressure  distribution  along  the  main  diagonal  of  the  computational  domain.  The 
computed  results  shown  are  for  the  finest,  i.e.  ,  21  x  21  grid.  As  is  evident,  the 
computed  results  are  in  excellent  agreement  with  the  exact  solution  for  all  values 
of  Reynolds  number. 

Figure  3.3(b)  shows  the  results  obtained  by  using  relatively  coarse  grids  for 
a  fixed  Reynolds  number  of  1.  It  can  be  seen  that  the  results  are  fairly  good  even 
at  rather  sparse  grids.  The  error  in  the  computed  pressure  distribution  is  higher 
at  points  close  to  the  boundary  of  the  domain  compared  to  the  errors  at  points  which 
are  in  the  core  of  the  domain. 

3.2  The  driven  cavity  problem 

The  problem  of  flow  in  a  square  cavity  with  one  moving  wall  is  sketched  in 
Fig.  3. A.  As  shown  in  this  figure,  three  sides  of  the  square  cavity  are  at  rest 
while  the  fourth,  the  top  wall,  moves  tangentially  with  a  velocity  u^.  Burgraff 
[4]  solved  this  problem  using  a  very  fine  (52  x  52)  grid.  Burgraff's  results  are 
taken  as  a  standard  of  comparison  for  the  results  obtained  using  the  proposed  equal 
order  method. 

Computational  details 

The  computational  domain  is  shown  in  Fig.  3.5.  The  domain  discretization  is 
also  shown  in  this  figure  for  a  13  x  13  grid.  The  grid  was  kept  finer  near  the 
walls  than  in  the  core  so  that,  at  higher  Reynolds  numbers,  the  thin  boundary  layers 
close  to  the  walls  could  be  resolved.  The  position  of  the  first  two  grid  lines 
parallel  to  each  wall  was  always  kept  fixed,  being  at  a  distance  of  .05  and  .1  from 
the  wall.  Grid  refinement  only  changed  the  spacing  6  between  the  interior  grid  lines 


Consequently,  for  a  N  x  N  grid, 


!-(.!)-(. 1)  .  -8 

Ws)  WT) 


(3.18) 


Results 

(i)  Figure  3.6(a)  shows  results  for  the  u  velocity  along  the  vertical  center- 
line  (x  =  L/2).  The  results  of  the  equal  order  method  are  for  a  21  x  21  (N=21) 
grid.  At  Reynolds  numbers  of  1  and  100,  the  computed  results  show  very  good  agree¬ 
ment  with  the  results  of  Burgraff;  in  particular,  the  difference  between  the  two 
cannot  be  resolved  on  the  graphical  scale  of  Fig.  3.6(a).  At  a  high  Reynolds 
number  of  400,  the  computed  results  are  infair  agreement  with  the  results  of  Burgraff 
except  close  to  the  region  where  the  gradient  of  velocity,  i.e.  3u/3y,  changes  sign. 
It  is  possible  that  the  21  x  21  grid  used  for  the  equal  order  method  is  not  enough  to 
resolve  the  thin  boundary  layers  close  to  the  walls  occurring  for  Re  =  400. 

Figure  3.6(b)  shows  how  the  solution  evolves  with  grid  refinement  by  the  equal 
order  method.  The  three  grid  sizes  shown  correspond  to  N  =  9 ,  13  and  21  with  cor¬ 
responding  6  (see  Fig.  3.5)  of  .2,  .1  and  .05. 

(ii)  For  two-dimensional  incompressible  flow,  the  stream  function  v:  is  defined 


such  that : 


from  which  it  follows 


u  ° 


v  * 


32^  .  92^ 

2  2 
3x*  3y* 


3ij> 

(3.19) 

3y 

3\1) 

3x 

(3.20) 

3u 

3v 

(3.21) 

*  3y  ' 

3x 

After  the  velocity  field  had  been  computed,  the  Poisson  equation  (3.21)  was  solved 
numerically  to  get  the  stream  function  at  all  grid  points. 
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The  streamline  plots  are  shown  in  Figs.  3.7  -  3.9  for  Reynolds  number  of  1, 

100  and  400  respectively.  The  results  computed  by  the  proposed  equal  order  method 
show  a  good  agreement  with  the  results  of  Burgraff.  The  separated  flows  at  the 
bottom  corners  are  correctly  picked  up  at  even  coarse  grids. 

3.3  Natural  convection  in  a  square  cavity 

The  problem  of  natural  convection  in  a  square  cavity,  shown  schematically  in 
Fig.  3.10,  is  used  as  another  test  problem.  The  hot  and  cold  vertical  walls  are 
both  isothermal  at  temperature  T^  and  T^,  respectively.  The  top  and  bottom  horizontal 
walls  are  adiabatic.  The  flow  is  assumed  to  be  incompressible  except  for  the  calcu¬ 
lation  of  driving  buoyancy  forces.  Also,  all  properties  of  the  fluid  are  assumed 
constant. 

Since  the  problem  does  not  have  an  exact  analytic  solution,  the  results  computed 
by  the  proposed  equal  order  method  will  be  compared  to  those  obtained  by  a  finite 
difference  technique  on  a  very  fine  (32  x  32)  grid.  Because  of  the  fine  grid  used, 
the  results  of  the  finite-difference  method  are  regarded  as  "exact".  The  domain  was 
discretized  into  a  uniform  N  x  N  grid  as  shown  in  Fig.  3.11. 

Results 

The  problem  is  governed  by  two  dimensionless  parameters,  namely,  the  Grashof 

number  and  the  Prandtl  number.  In  the  presented  computations,  the  value  of  Prandtl 

3  4 

number  was  kept  fixed  at  1,  while  the  Grashof  number  was  assigned  values  10  ,  10 
and  10^ . 

(i)  Figures  3.12-3.14  dhow  the  variation  of  vertical  velocity  and  temperature 
along  the  horizontal  mid  plane  (y/L  •  .5)  for  Grashof  numbers  of  10  ,  104'  and  10  . 
Even  though  the  results  of  the  proposed  method  have  been  obtained  on  a  coarser  grid, 
their  accuracy  is  as  good  as  that  given  by  the  fine-grid  finite-difference  solution. 

(ii)  The  average  Nusselt  number  on  the  hot  wall  is  defined  as 


Fig.  3.10  Natural  convection  in  a  cavity 


Fig.  3.11  Domain  discretization 
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h  L 

Nu  -  (3.22) 

a  v  k 

where,  h  is  the  average  heat  transfer  coefficient  defined  from 

av 

qav 

hav  "  (Th-Tc)  (3.23) 

qav,  the  average  heat  flux  being  given  by, 

q  -  f  /  q(y)dy  (3.24) 

Hav  L  0 

Table  3-3  shows  the  results  for  the  average  Nusselt  number  on  the  hot  wall 
computed  using  the  equal  order  method.  The  results  have  been  cross  tabulated 
against  the  grid  size  and  Grashof  number.  For  extrapolating  the  computed  results 
to  those  corresponding  to  6  -  0  (infinitely  fine  grid),  a  Richardson  extrapolation 
has  been  applied.  Thus 


Nuav(<5)  -  Nugv(0)  +  a  6P 


(3.25) 


Table  3'3 

Average  Nusselt  number  in  the  duct,  convergence  study 


Grid 

5 

Nu  ..  for  Gr  » 
aV 

103 

10- 

105 

11  x  11 

.1 

1.1084 

2.2112 

4.2611 

15  x  15 

.07143 

1.1125 

2.2302 

4.4069 

19  x  19 

.05556 

1.1146 

2.2408 

4.4719 

Extrapolated  value 
as  6  ■* 

1.1202 

2.2793 

4.5909 

Value  obtained  by 

F  0  method;  32  *  32 
Grid 

1.1203 

2.2825 

4.7755 
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where  a  and  p  are  some  constants,  and  Nu  (0)  is  the  value  of  Nu  for  6=0.  If 

av  av 

Nu^v  is  computed  for  three  values  of  6,  Equation  (3.25)  can  be  inverted  to  obtain 

Nu  (0),  a  and  p.  Table  3-3  includes  the  value  of  Nu  (0). 
av  r  av 

It  is  clear  from  Table  3-3  that  the  computed  values  of  average  Nusselt  number 
are  in  very  good  agreement  with  those  obtained  by  the  fine-grid  finite-difference 


method. 
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where  the  diffusion  term  in  the  main-stream  direction  z  has  been  omitted  in 
conformity  with  the  boundary-layer  (or  parabolic-flow)  assumption.  Once  again, 
the  general  dependent  variable  <p  can  stand  for  enthalpy,  concentration,  turbulence 
parameters,  and  the  velocity  components,  u,  v,  and  w.  When  <j>  equals  u  or  v,  the 
source  terms  S  would  contain  the  pressure  gradient  -  dp/dx  or  -  3p/3y;  these  will 
be  treated  by  the  procedure  described  in  Section  2.  When  <p  stands  for  the  main¬ 
stream  velocity  w,  the  source  term  would  include  -  dp/dz.  However,  this  will  be 
replaced  by  -  dp/dz,  where  p  represents  a  cross-sectional  mean  pressure.  To 
regard  the  velocity  w  to  be  "driven"  by  a  mean  pressure  p  and  not  affected  by  the 
local  variation  of  pressure  p  constitutes  the  basis  of  the  fully  parabolic  treat¬ 
ment.  If  the  velocity  w  is  allowed  to  be  influenced  by  the  local  pressure  p,  the 
resulting  procedure  is  called  partially  parabolic.  Construction  of  a  partially 
parabolic  procedure  represents  a  later  task  in  the  current  research  project. 

4.2  Domain  discretization 

The  advantage  of  a  marching  or  parabolic  procedure  is  that,  although  the 
flow  domain  is  three-dimensional,  the  entire  duct  need  not  be  considered  at  once. 

At  any  given  station,  the  computational  problem  is  to  obtain,  from  the  known 
values  of  the  variables  on  an  upstream  plane  (see  Fig.  4.1),  the  unknown  values 
of  the  variables  on  the  next  downstream  plane.  Successive  repetition  of  this  basic 
operation  is  used  to  cover  the  total  length  of  the  duct.  Restriction  of  the  basic 
computational  module  to  the  region  between  two  planes  implies  that  computer  storage 
is  needed  for  the  (j)  values  only  on  the  two  planes  and  not  throughout  the  entire 
duct. 

The  domain  is  discretized  by  forming  triangular  elements  over  the  upstream 
and  downstream  planes.  Corresponding  typical  elements  are  shown  in  Fig.  4.2.  The 
shape  of  the  two  planes  need  not  be  rectangular;  any  arbitrary  shape  of  the  duct 
cross  section  can  be  adequately  represenred  by  triangular  elements.  Further,  the 
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shape  and  size  of  the  downstream  plane  neec  not  match  the  upstream  plane. 

Indeed,  an  ability  to  handle  varying  cross  sections  of  arbitrary  shape  is  one  of 
the  objectives  of  the  proposed  calculation  method.  Only  when  the  cross  section 
changes  too  rapidly  from  one  station  to  the  next,  does  the  proposed  fully  parabolic 
procedure  become  inapplicable;  then  the  problem  should  be  solved  by  a  partially 
parabolic  or  even  a  fully  elliptic  procedure. 


Fig.  4.2  Typical  elements  on  the  upstream  and  downstream  planes 
The  triangulation  on  the  two  planes  must,  however,  have  the  same  number  of 
triangles  with  a  one-to-one  correspondence  between  the  elements.  Thus,  the  tri¬ 
angulation  of  the  downstream  plane  can  be  regarded  as  a  stretched,  compressed,  or 
distorted  version  of  the  triangulation  on  the  upstream  plane.  The  element  PKL 
shown  in  Fig.  4.2  corresponds  to  the  element  PU-KU-LU  on  the  upstream  plane:  the 
size  and  orientation  of  these  elements  need  not  be  the  same.  The  dashed  lines 
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show  how  each  triangle  is  subdivided  for  the  construction  of  the  control  volumes. 
The  upstream  and  downstream  faces  of  a  typical  control  volume  resemble  the  shaded 
region  shown  in  Fig.  2.1.  The  lateral  faces  of  the  control  volume  are  formed  by 
quadrilaterals  such  as  bu-b-d-du  formed  by  the  points  shown  in  Fig.  4.2.  The 
control  volume  has  thus  the  shape  of  a  prism  with  the  shaded  polygonal  region  in 
Fig.  2.1  as  its  base. 

4.3  Discretization  equation 

The  discretization  form  of  Eq.  (4.1)  is  obtained  by  integrating  the  equation 
over  a  typical  control  volume.  For  the  convection  and  diffusion  fluxes  in  the  x 
and  y  directions  (across  the  lateral  faces  of  the  control  volume),  the  treatment 
described  in  Section  2.3  including  the  exponential  shape  function  is  employed. 

The  z-direction  convection  across  the  upstream  and  downstream  faces  of  the  control 
volume  is  obtained  by  assuming  that  the  values  <f>pu  and  <j>p  respectively  prevail  over 
these  faces.  Since  it  is  possible  that  the  lateral  faces  are  not  exactly  parallel 
to  the  z  axis,  there  may  also  be  a  z-direction  convection  across  the  lateral  faces. 
This  is  calculated  in  terms  of  the  4)  value  obtained  at  the  control  volume  faces 
from  the  shape  function  (given  by  Eq.  (2.2)). 

In  expressing  the  fluxes  across  the  lateral  faces,  one  is  faced  with  the 
choice  of  using  the  (known)  values  of  <p  at  the  upstream  plane,  the  (unknown)  values 
of  4>  at  the  downstream  plane,  or  some  combination  of  these.  Thus,  the  choice  is 
between  the  explicit,  fully  implicit,  or  Crank-Nicolson  formulations,  or  any  other 
variants  of  these.  As  a  result  of  the  considerations  expressed  in  [3],  the  fully 
implicit  procedure  is  adapted  here.  Thus,  the  known  4>  values  on  the  upstream 
plane  are  used  only  in  the  calculation  of  the  convection  through  the  upstream  face 
of  the  control  volume.  All  other  influences  in  the  discretization  equation  are 
expressed  in  terms  of  the  unknown  <j>  values  on  the  downstream  plane. 
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The  final  discretization  equation  can  be  written  as 

aP  *P  *  Z  anb  ^nb  +  aPU  ^PU  +b 


(4.2) 


where  apy  is  the  mass  flow  rate  through  the  upstream  face  of  the  control  volume. 
With  the  values  known  from  the  upstream  plane,  Eq.  (4.2)  has  the  same  form 
as  Eq.  (2.3);  therefore,  the  solution  scheme  outlined  in  Section  2.4  is  appli¬ 
cable  to  Eq.  (4.2)  as  well. 

4.4  Calculation  of  the  overall  pressure  gradient 

In  a  duct  flow,  the  pressure  gradient  dp/dz  driving  the  z-direction  velocity 
w  is  not  known  beforehand.  It  must  adjust  itself  so  that  a  given  mass  flow  rate 
M  can  be  maintained  in  the  duct  despite  the  resistance  offered  by  the  duct  walls. 
Thus,  the  satisfaction  of  the  overall  mass  conservation  is  a  constraint  that  per¬ 
mits  the  determination  of  dp/dz. 

The  discretization  equation  for  the  mainstream  velocity  w  can  be  written  along 
the  lines  of  Eq.  (4.2).  It  is 

‘p,  “p  -  £  anb  "nb  +  “TO  "pu  +  b 


-  (AV)  (dp/dz) 


(4.3) 


where  the  pressure-gradient  term  is  written  separately;  its  multiplier  AV  is  the 
volume  of  the  control  volume.  The  determination  of  dp/dz  is  based  on  the  procedure 


suggested  by  Raithby  and  Schneider  [5].  The  main  idea  of  the  procedure  is  that, 
for  a  given  set  of  coefficients,  all  the  w  values  on  the  downstream  plane  are 
linearly  dependent  on  dp/dz.  Thus 

wp  -  ctp  -  6p  (dp/dz)  (4.4) 

where  Op  and  Bp  are  unknown  constants  associated  with  each  node  P.  Equation  (4.4) 
can  be  substituted  into  Eq.  (4.3).  Since  the  resulting  relation  must  hold  for  any 
value  of  dp/dz,  two  separate  equations  can  be  derived:  one  composed  of  terms  that 
do  not  contain  dp/dz  and  the  other  composed  of  the  coefficients  of  dp/dz.  Thus, 


*P  “P  ^  nb  nb  '  “PU  "PU 


and  — * 

*p  6p  ■  Ea„b  &>  * 4V  (4-«  ^ 

These  equations  have  the  same  general  form  of  Eq.  (4.2)  and  can  be  solved  in  a 

similar  fashion  to  obtain  ap  and  8p  for  every  node. 

The  total  mass  flow  rate  M  through  the  duct  is  given  by 

M  =  £  pp  A?  wp  (4.7) 

where  pp  and  Ap  denote  the  density  and  area  associated  with  the  downstream  face  of 
the  control  volume  and  the  summation  is  taken  over  the  entire  downstream  plane. 

Combination  of  Eqs.  (4.4)  and  (4.7)  gives 

d?/dz  -  (ft  -  Epp  Ap  ap)/(Epp  Ap  6p)  (4.8) 

4.5  The  overall  solution  procedure 

The  complete  solution  of  a  three-dimensional  duct  flow  is  obtained  by  repeating 
the  solution  for  one  forward  step  in  the  z  direction.  For  the  first  forward  step, 
the  values  of  <p  on  the  upstream  or  the  inlet  plane  are  given  as  a  part  of  the  problem 
specification.  For  subsequent  forward  steps,  the  <J>  values  obtained  on  the  downstream 
plane  of  the  previous  step  become  available  as  the  upstream-plane  values  for  the 
current  step.  With  this  general  framework,  it  is  now  sufficient  to  describe  the 
details  of  the  calculation  for  one  forward  step.  The  various  steps  in  the  calculation 
sequence  are  outlined  here. 

(1)  Start  with  the  initial  guess  for  the  <f>  values  for  the  downstream  plane. 

Usually,  the  known  <j>  values  on  the  upstream  plane  serve  as  satisfactory  guesses. 

(2)  Calculate  the  discretization  coefficients  for  the  three  momentum  equations. 

(3)  Using  these  coefficients,  solve  Eqs.  (4.5)  and  (4.6)  to  obtain  the  fields 
of  a  and  6. 

(4)  Determine  dp/dz  from  Eq.  (4.8)  and  hence  obtain  the  values  of  w  from 


Eq.  (4.4). 
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(5)  Now  calculate  u,  v,  du,  and  dV  from  Eqs.  (2.7),  (2.8),  and  (2.11)  appropri¬ 
ately  written  for  the  three-dimensional  control  volume. 

(6)  Hence  solve  the  three-dimensional  counterpart  of  Eq.  (2.14)  to  obtain  the 
pressure  field. 

(7)  Considering  this  as  p*,  solve  the  cross-stream  momentum  equations  to  get  u* 
and  v*  . 

(8)  Now  set  up  the  pressure-correction  equation  (2.16)  and  solve  it.  Use  the 

resulting  p'  values  to  correct  u*,  v*,  u*,  and  v*  . 

(9)  Solve  Eq.  (4.2)  for  other  relevant  variables  such  as  temperature,  concen¬ 
tration,  turbulence  parameters,  etc. 

(10)  Return  to  step  2  and  repeat  until  convergence. 

(11)  Treat  the  downstream  0  values  as  the  upstream  values  for  the  next  forward 
step  and  return  to  step  1  to  begin  the  calculation  sequence  for  the  next  Az. 
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5.  Application  of  the  method  to  some  three-dimensional  duct  flows 

In  this  section,  results  obtained  by  applying  the  proposed  method  to  sense 
three-dimensional  ducts  flows  are  presented.  The  results  are  divided  into  two 
sub-scetions:  In  Section  S.l,  ducts  of  uniform  cross-sectional  shape  and  area 
are  considered;  in  Section  5.2,  ducts  of  variable  cross  section  are  treated. 

5.1  Flow  through  ducts  of  uniform  cross  section 

Here,  two  duct  geometries  are  considered.  The  first  example  involves  dev¬ 
eloping  flow  in  a  square  duct.  The  second  exanple  considers  the  flow  over  rod 
bundles.  This  problem  demonstrates  the  capability  of  the  proposed  method  to 
handle  flow  through  ducts  of  complex  (irregular)  cross-sectional  shapes. 

5.1-1  Developing  flow  in  a  square  duct 

The  flow  situation  considered  is  shown  schematically  in  Fig.  5.1.  A  uniform 
flow  enters  a  square  duct  at  the  entrance  plane  z  *  0.  As  the  flow  moves  along 
the  duct,  boundary  layers  on  the  duct  walls  grow  continuously.  After  a  certain 
distance,  they  merge  into  each  other  and  the  flow  tends  to  attain  a  fully  developed 
character.  Such  a  fully  developed  flow  is  characterized  by  zero  cross-stream 
velocities,  i.e.  ,  u  ■  v  -  0,  and  an  axial  velocity  w  which  is  independent  of  the 
axial  location  z.  Further,  in  fully  developed  flow,  the  axial  pressure  gradient 
dp/dz  becomes  independent  of  z  and  hence  p  decreases  linearly  with  z. 
Computational  details 

Because  of  symmetry,  only  one  fourth  of  the  duct  is  used  for  computational 
purposes.  Figure  5.2  shows  the  cross  section  of  one  fourth  of  the  duct  and  its 
discretization  into  triangular  elements.  As  can  be  seen,  a  uniform  grid  having 
N  x  N  grid  points  is  used.  For  a  N  x  N  grid,  the  spacing  between  adjacent  grid 
points  is  given  by  6  where 


Results 


(1)  The  dependence  of  the  computed  results  on  the  grid  parameters  Az  and 

6  Is  shown  by  reference  to  fRe,  where  f  is  the  friction  factor  defined  as 

D 


(-  - 

'  dz'  1 . 


and  Re  is  the  Reynolds  number  given  by 


Re 


dz' 

pwD 

U 


2^w 


(5.2) 


(5.3) 


Figure  5.3  shows  the  computed  variation  of  fRe  with  dimensionless  axial 
distance.  Results  obtained  by  using  different  step  sizes  Az  and  different  number 
of  grid  points  in  the  cross-stream  plane  are  shown.  The  finest  grid  used  for  the 
results  displayed  in  Fig.  5.3  corresponds  to  a  dimensionless  step  size  of  .0005 
and  a  17  x  17  grid  in  the  cross-stream  plane.  All  the  remaining  results  to  be 
presented  are  based  on  this  grid. 

(ii)  The  variation  of  mean  pressure  p  with  axial  distance  z  is  shown  in 
Fig.  5.4(a).  Also  shown  are  the  experimental  results  of  Beavers,  et  al  [6], 

The  agreement  between  the  computed  and  experimental  results  is  very  good. 

Figure  5.4(b)  shows  the  variation  of  axial  velocity  at  the  centerline  of  the 
duct.  The  experimental  results  of  Goldstein  and  Kreid  [7]  have  also  been  plotted 
for  comparison.  Once  again,  the  agreement  between  the  computed  and  experimental 
results  is  very  good. 

5.1-2  Developing  flow  over  rod  bundles 

The  second  test  problem  involves  computation  of  developing  flow  over  rod 
bundles  which  are  arranged  in  a  regular  equilateral  triangular  array.  The  situ¬ 
ation  is  shown  schematically  in  Fig.  5.5.  Such  a  flow  configuration  arises  in 
nuclear  reactors,  compact  heat  exchangers  and  many  other  engineering  applications 
Computational  details 


Because  of  symmetry,  the  calculation  can  be  confined  to  the  region  shown  in 
Fig.  5.6.  The  corresponding  triangulation  is  also  shown.  A  17  x  17giid  was  used 


Fig. 
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Figure  5.7(a)  shows  the  variation  of  pressure  p  with  axial  distance 
for  different  values  of  the  geometrical  parameter  s/rQ.  The  development  of 
centerline  velocity  is  shown  in  Fig.  5.7(b).  The  fully  developed  values  of  the 
centerline  w  are  in  good  agreement  with  thos  obtained  by  Sparrow  and  Loeffler  [8). 
5.2  Flow  through  ducts  of  varying  cross-sectional  shape  and  area 

The  most  important  feature  of  the  proposed  method  is  its  capability  to 
compute  the  flow  through  ducts  of  varying  cross-sectional  shape  and  area.  Since 
no  experimental  or  analytical  results  are  available  for  such  ducts,  a  procedure  is 
devised  here  to  "fabricate"  exact  solutions  for  developing  flow  in  ducts  of  varying 
cross  sections.  These  exact  solutions  are  used  to  check  the  accuracy  of  the  numeri¬ 
cal  solution. 

5.2-1  Formulation  of  the  exact  solutions 

The  central  idea  of  the  method  can  be  summarized  as  follows: 

(i)  Propose  a  velocity  field  u  which  satisfies  the  continuity  equation 

6XaC  t 

and  certain  boundary  conditions.  Also  propose  a  pressure  field. 

U  V  V 

(il)  Compute  the  source  terms  in  the  momentum  equations,  S  ,  S  ,  and  S  ,  such 

-4 

that  the  velocity  field  u  satisfies  the  momentum  equations. 

exact 

If  the  resulting  expressions  for  the  source  terms  are  considered  as  given 

”4 

"body  forces"  in  the  problem  specification,  then  uexact  is  indeed  the  exact  solution 
of  the  flow  equations  (along  with  the  assumed  pressure  field) . 

Case  A:  Square  diffuser  with  parabolic  velocity  profile 

For  a  square-sectioned  diffuser,  shown  in  Fig.  5.8,  the  half-side  D  linearly 
increases  at  a  rate  of  dD/dz  ■  a.  The  exact  solution  for  the  velocity  field  is 
given  by 

,  -  9  H  (P2-y2)(P2-x2)  (5.4) 

16  o  ^5 
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where  M  is  the  mass  flow  rate  through  the  duct. 


Case  B:  Elliptic  ducts 
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Fig.  5.9  Cross  section  of  an  elliptic  duct 

The  cross  section  of  an  elliptic  duct  is  shewn  in  Fig.  5.9.  The  semi-major 


axis  D  and  the  semi-minor  axis  D  are  both  functions  of  the  axial  coordinate  z. 

y  x 


Let  the  variation  of  D  and  D  be  given  by: 

x  y 


°0  *  ~<z 


(5.13) 


D  D  ■  constant  *  D- 
y  x  0 


(5.14) 


where  Y  and  DQ  are  two  arbitrary  constants.  Thus,  at  z  =  0 ,  the  cross-section  of 
the  duct  is  a  circle  of  radius  Dq.  The  exact  solution  for  the  velocity  field  can 
be  expressed  as 
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provided  that  the  source  terms  in  the  momentum  equations  are  given  by 


(5.21) 


S 
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_8  x§_ 

*2  p  D1? 


y[(Dx2-x2)+(D2-y2)] 


(5.22) 


Here,  the  coefficient  6  is  defined  as 


(5.24) 


In  the  use  of  the  calculation  procedure  for  these  test  problems,  computations 
were  performed  both  with  and  without  the  momentum  source  terms.  The  results  of 
the  former  can  be  comared  with  the  exact  solutions,  while  the  latter  computation 
shows  how  the  flow  will  behave  in  reality  without  the  artificial  "body  forces". 
5.2-2  Flow  through  a  sauare  diffuser 


With  reference  to  Fig.  5.8,  the  geometry  of  the  square  diffuser  is  defined 
in  terms  of: 


Thus, 


2D  -  D.  +  ctz 
i  n 


(5.25) 


(5.26) 


D  -  DQ  +  yz 


(5.27) 


where 


D0  "  Din/2 


Y  *  a/2 


(5.28) 

(5.29) 


The  computations  were  performed  by  using  a  17  x  17  grid  similar  in  pattern  to 
the  one  in  Fig.  5.2. 

Figure  5.10  shows  the  variation  of  axial  pressure  gradient  dP/dZ  with  Z, 


where  the  dimensionless  distance  Z  is  defined  by 

Z  -  z/ (Re .  D,  )  (5.30) 

m  m 

Figure  5.10(a)  shows  the  results  with  the  source  terms,  while  Fig.  5.10(b) 
displays  the  implications  of  the  usual  momentum  equations.  The  dashed  lines  in 


Fig.  5. 10(a)  represent  the  exact  solution;  the  agreement  can  be  seen  to  be  very 
good.  For  the  flow  without  the  extra  source  terms,  separation  is  encountered  for 
large  angles  of  the  diffuser. 

The  variation  of  centerline  velocity  with  the  axial  distance  Z  is  shown  in 
Fig.  5.11.  Once  again,  good  agreement  with  the  exact  solution  is  seen  in 
Fig.  5.11(a),  while  the  flow  separation  in  absence  of  the  special  source  terms  is 
evident  in  Fig.  5.11(b).  The  development  of  the  axial  velocity  profiles  for  the 
case  of  zero  source  terms  is  shown  in  Fig.  5.12.  The  profiles  at  flow  separation 
can  be  seen  to  exhibit  an  inflexion  point  at  the  wall. 

5.2-3  Flow  through  an  elliptic  duct 

For  the  flow  through  an  elliptic  duct,  the  following  geometrical  relations 
are  used. 


20  “  D;„  +  az 
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°y  ■  D0/0x 

(5. 

Thus,  at  the  inlet,  the  cross  section  of  the  duct  is  a  circle  of  diameter 

(radius  Dn  »  D.  /2) . 
u  in 

Because  of  symmetry,  only  one  quarter  of  the  duct  needs  to  be  considered.  A 


cross  section  of  the  quarter  of  the  duct,  and  its  discretization  into  triangular 
elements  is  shown  in  Fig.  5.13.  A  17  x  17  grid  is  used. 


The  variation  of  dimensionless  axial  pressure  gradient  dP/dZ  with  the  axial 


distance  Z  is  shown  in  Fig.  5.14  for  a  range  of  values  aRe.  .  The  results  for  the 

in 

flow  with  source  terms  are  displayed  in  Fig.  5.14(a)  and  the  results  without  the 

source  terms  are  shown  in  Fig.  5.14(b).  It  is  evident  from  Fig.  5.14(a)  that 

the  computed  results  are  in  good  agreement  with  the  exact  solution  even  for  large 

values  of  aRe.  . 

in 

The  variation  of  axial  velocity  (w/w.  )  on  the  semi-major  and  semi-minor  axes 

in 

of  the  ellipse  for  the  flow  with  source  terms  is  shown  in  Fig.  5.15.  Also  shown 
is  the  exact  solution.  The  agreement  between  the  computed  and  exact  velocity 
profiles  can  be  seen  to  be  very  good.  The  velocity  profiles  for  the  flow  without 
source  terms  are  shown  in  Fig.  5.16. 

6.  Concluding  remarks 

This  report  has  described  the  development  of  a  calculation  procedure  for 
three-dimensional  flow  in  ducts  of  varying  cross  sections.  The  procedure  is  dev¬ 
eloped  in  two  stages.  First,  a  calculation  scheme  for  two-dimensional  flows  is 
described.  Then,  it  is  used  as  the  solution  method  over  a  cross  section  of  a  duct 
in  a  marching  procedure.  Examples  of  both  two-dimensional  flows  and  three-dimen¬ 
sional  duct  flows  are  provided  to  demonstrate  the  capabilities  of  the  proposed 
procedure. 

Further  work  in  this  research  project  consists  of  additional  testing  of  the 
method  and  the  development  of  pratially  parabolic  and  fully  elliptic  procedures 
for  more  complex  duct  flows. 
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